About me

• I earned my PhD in Evolutionary Anthropology during the 433rd lunation of the third millennium, when Mars was at an orbital longitude of ~311 degrees in May 2022

• To get from Zurich to here, I travelled 807 light-microseconds 242 kilometers

• The units we use to communicate matter!

The way we report statistics is often confusing (and misleading)

The way we report statistics is often confusing (and misleading)

In contrast, descriptive statistics and data viz are easy to understand

‘Monolingual participants have a mean cognitive flexibility score of 51.27 (SD = 18.13), while multilingual participants have a mean score of 61.77 (SD = 18.78).’

• Simple comparisons with scientifically-meaningful units.
• What if we could do the same for statistical models?

Predictions, not parameters

• The central idea: rather than directly interpret parameters, use the model to make predictions
  • contrasts between different predictions consitute marginal effects
• The marginaleffects package by Vincent Arel-Bundock makes these calculations easy for a wide range of models in both R and Python.

In simple models, a parameter might be all we need

Call:
lm(formula = cognitive_flexibility ~ multilingual)

Residuals:
    Min      1Q  Median      3Q     Max 
-52.387 -11.913   1.581  13.388  36.985 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)    51.267      1.665  30.795  < 2e-16 ***
multilingual   10.508      2.336   4.499 1.05e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 18.46 on 248 degrees of freedom
Multiple R-squared:  0.07545,   Adjusted R-squared:  0.07172 
F-statistic: 20.24 on 1 and 248 DF,  p-value: 1.051e-05

A regression is a machine that makes predictions

\[ \text{Cognitive flexibility} = \alpha + \beta \times \text{Multilingual} + \epsilon \] \[ \epsilon \sim \text{Normal}(0, \sigma^2) \]

• The parts of the machine are the parameters (\(\alpha, \beta, \sigma\)):
  • To make a prediction for a monolingual child, just take the intercept (\(\alpha\)).
  • To make a prediction for a multilingual child, add the coefficient for multilingual (\(\beta\)) to the intercept.
• Linear regression is such a simple machine that it offers a 1:1 mapping between the parameters and the estimand (average difference between multilingual and monolingual children).

Interactions

Interactions imply heterogeneity

Interactions imply heterogeneity

Call:
lm(formula = cognitive_flexibility ~ multilingual * age_c)

Residuals:
    Min      1Q  Median      3Q     Max 
-45.937  -9.753   1.457  12.291  38.200 

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)          51.123      1.580  32.353  < 2e-16 ***
multilingual         10.969      2.217   4.948 1.39e-06 ***
age_c                 2.316      1.079   2.145   0.0329 *  
multilingual:age_c    2.949      1.506   1.959   0.0512 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 17.51 on 246 degrees of freedom
Multiple R-squared:  0.1752,    Adjusted R-squared:  0.1652 
F-statistic: 17.42 on 3 and 246 DF,  p-value: 2.738e-10

• Which coefficient corresponds to the estimand (average marginal effect)?
  • None! The AME depends on the distribution of the moderator (age):
    • The magnitude of the multilingual effect depends on age–now it matters which children we are averaging over.

Meaning of interaction term is highly contextual

On the previous slide, the interaction term was mean-centered. What if we use the raw age variable instead?

Call:
lm(formula = cognitive_flexibility ~ multilingual * age)

Residuals:
    Min      1Q  Median      3Q     Max 
-45.937  -9.753   1.457  12.291  38.200 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)   
(Intercept)        29.321     10.350   2.833   0.0050 **
multilingual      -16.799     14.351  -1.171   0.2429   
age                 2.316      1.079   2.145   0.0329 * 
multilingual:age    2.949      1.506   1.959   0.0512 . 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 17.51 on 246 degrees of freedom
Multiple R-squared:  0.1752,    Adjusted R-squared:  0.1652 
F-statistic: 17.42 on 3 and 246 DF,  p-value: 2.738e-10

Getting started with marginaleffects: predictions

library(marginaleffects)

uncentered_model <- lm(cognitive_flexibility ~ multilingual * age, data = d)

preds <- predictions(uncentered_model)

preds |> 
    as.data.frame() |> 
    dplyr::select(rowid, multilingual, age, estimate, conf.low, conf.high) |> 
    knitr::kable(digits = 2)

Getting started with marginaleffects: predictions

library(marginaleffects)

centered_model <- lm(cognitive_flexibility ~ multilingual * age_c, data = d)

preds <- predictions(centered_model)

preds |> 
    as.data.frame() |> 
    dplyr::select(rowid, multilingual, age_c, estimate, conf.low, conf.high) |> 
    knitr::kable(digits = 2)

From predictions to comparisons

comparisons <- comparisons(uncentered_model, variables = "multilingual")

comparisons |> 
    as.data.frame() |> 
    dplyr::select(rowid, contrast, age, estimate, conf.low, conf.high) |> 
    knitr::kable(digits = 2)

Calcualting average marginal effects

AME <- avg_comparisons(uncentered_model)

AME |> 
    as.data.frame() |> 
    dplyr::select(term, contrast, estimate, conf.low, conf.high) |> 
    knitr::kable(digits = 2)

marginaleffects can also make plots

plot_predictions(uncentered_model, condition = list("multilingual"))

marginaleffects can also make plots

plot_predictions(uncentered_model, condition = list("age"))

marginaleffects can also make plots

plot_predictions(uncentered_model, condition = list("age", "multilingual"))

marginaleffects can also make plots

Because the plots are made with ggplot2, we can customize them as we like.

plot_predictions(uncentered_model, condition = list("age", "multilingual")) +
    theme(legend.title = element_blank()) +
    scale_color_manual(values = c("monolingual" = "royalblue", 
                                "multilingual" = "darkorange")) +
    scale_fill_manual(values = c("monolingual" = "royalblue", 
                                "multilingual" = "darkorange")) +
    labs(x = "Age", y = "Cognitive flexibility score")

Generalized linear models

Logistic regression

• Hypothetical dataset: Touchscreen experiment with Chimpanzees
  • After viewing a stimulus, the chimpanzee touches the screen on either the ‘Agent’ or the ‘Patient’ in a scence.
  • Estimand: what is the average difference in the probability of touching the ‘Agent’ side between two conditions (‘experimental’ vs ‘control’)?
    • \(\frac{1}{N} \sum_{i=1}^N [P(Y_i = \text{Agent} | \text{experimental}) - P(Y_i = \text{Agent} | \text{control})]\)

Logistic regression

\[ \text{Agent} \sim \text{Bernoulli}(\pi) \] \[ \pi = \text{logit}^{-1}(\alpha + \beta \times \text{experimental}) \]

Call:
glm(formula = agent ~ condition, family = binomial(link = "logit"), 
    data = d)

Coefficients:
                      Estimate Std. Error z value Pr(>|z|)  
(Intercept)             0.5465     0.3789   1.443   0.1491  
conditionexperimental   1.3253     0.6573   2.016   0.0438 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 67.48  on 59  degrees of freedom
Residual deviance: 62.99  on 58  degrees of freedom
AIC: 66.99

Number of Fisher Scoring iterations: 4

Logistic regression

• Even with just a single predictor, we cannot get our estimate from looking at model coefficients.

• Due to the link function, the effect of the experimental condition depends on the intercept.

• ‘Mechanistic’ rather than conceptual interaction.

  • To get the AME we need to calculate:
    • \(\frac{1}{1 + e^{-(\alpha + \beta \times \text{experimental})}} - \frac{1}{1 + e^{-\alpha}}\)

Logistic regression

The marginal effect on the probability scale for a log odds ratio of 1.33 depends strongly on the intercept!

marginaleffects to the rescue

AME <- avg_comparisons(logistic_model, variables = "condition")

AME |> 
    as.data.frame() |> 
    dplyr::select(term, contrast, estimate, conf.low, conf.high) |> 
    print(digits = 2)

       term                           contrast estimate conf.low conf.high
1 condition mean(experimental) - mean(control)     0.23    0.022      0.44

plot_predictions(logistic_model, condition = list("condition")) + ylab("Pr(Agent)")

• NOTE: I should write out the estimate in line with the original goal to make the outputs easier to report and interpret.

But how does it work?

• For frequentist models, marginaleffects by default uses the delta method (Taylor series expansion)

  • Various bootstrap and weighting methods supported

• For Bayesian models (e.g., fit with brms), marginaleffects, uses full posterior draws

• More info: https://marginaleffects.com/chapters/uncertainty.html

Review thus far

Reporting model results using marginaleffects:

• Makes your model results easier to communicate
• Avoid misinterpreting the magnitude of effects
• Avoid misinterpreting the direction of effects

Non-linear effects

Example: relationship between population density and structural diversity in language (measured by normalized entropy)

Non-linear effects

Family: gaussian 
Link function: identity 

Formula:
structural_diversity ~ s(log_density, k = 8)

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 0.676467   0.009797   69.05   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Approximate significance of smooth terms:
                 edf Ref.df     F p-value    
s(log_density) 6.221  6.798 49.08  <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

R-sq.(adj) =  0.623   Deviance explained = 63.5%
GCV = 0.019914  Scale est. = 0.019195  n = 200

Non-linear effects

slopes_df <- slopes(gam_model) |> 
    as.data.frame() |> 
    dplyr::select(estimate, conf.low, conf.high, log_density, structural_diversity)

print(head(slopes_df), digits = 2)

  estimate conf.low conf.high log_density structural_diversity
1   0.1590    0.110     0.208        0.88                 0.45
2   0.0023   -0.059     0.064        5.88                 0.89
3  -0.0065   -0.059     0.046        2.09                 0.68
4   0.0392   -0.020     0.098        6.83                 0.83
5   0.0912   -0.015     0.197        7.40                 0.96
6   0.0240   -0.100     0.148       -1.54                 0.37

Non-linear effects: average slopes

\(\mathbb{E}\left[\frac{\partial \text{structural diversity}}{\partial \text{ log pop density}}\right]\)

avg_slopes(gam_model) |> 
    as.data.frame() |> 
    dplyr::select(term, estimate, conf.low, conf.high) |> 
    knitr::kable(digits = 2)

The bigger picture

• Directly interpreting statistical models is a full of traps that can mislead us. Working with predictions and marginal effects allows us to avoid many of these traps.

• marginaleffects can be used with a huge array of models:

  • linear models, generalized linear models (e.g., lm, glm)
  • generalized additive models (e.g., mgcv, gamm4)
  • multilevel models (e.g., lme4, glmmTMB)
  • Bayesian models (e.g., brms)
  • machine learning models (e.g., tidymodels, mlr3)

Can marginal effects improve the reliability of science?

• On the one hand, we can avoid falling into many of the traps of interpreting model parameters.

• On the other hand, this way of doing things brings the estimand into focus. We are being more explicit about the quantity we are trying to estimate–which is not tied to any particular model.

Red Card Dataset

2,053 soccer players in the first male divisions of England, Germany, France, and Spain in the 2012–2013 season.

Outcome: number of red cards recieved for 146,028 player-referee dyads.

Covariates: player and referee traits, match statistics

Many analysts, many estimates

“Are soccer players with dark skin tone more likely than those with light skin tone to receive red cards from referees?”

Many analysts, many estimates

“Are soccer players with dark skin tone more likely than those with light skin tone to receive red cards from referees?”

Who watches the watchmen? Reanalysis of the meta-analysis

Many analysts, many estimands

“Thus, we argue that the results obtained in the CSI showed such a large variation because the 29 teams pursued (at least) four different research questions and therefore used different research designs.”

Controlled direct effect: \(E[Y(\text{dark}) - Y(\text{light}) | M = m]\); 486 parametric models

Causal (not statistical) assumptions needed to inform research design

Breakout session (10-15 minutes)

Report-back/Q&A session

• How did it go? What friction did you encounter?

Statistics can feel like a thankless list of things that all have to be right or else nothing is right. - Not totally wrong, but leads to anxiety

• Calls to improve statistical practice in science imply time commitments:
  • More time for training
  • More time for coding tasks (recall the multiverse)

Academic time allocation

Brauer, K. (2021). ” I’ll Finish It This Week” And Other Lies. arXiv preprint arXiv:2103.16574.

Harsh realities

• Time allocation is a zero-sum game
  • Time spent on doing better statistics is time that could be spent on the many other aspects of science and professional development
• If raise the bar on statistical pracitce we need workflows that are reproducible, understandable, and which reduce cognitive burden.

The tangled web of data analysis

• What happens when Reviewer #2 asks for a change to the data/model?

targets

Package for reproducible research pipelines by Will Landau:

• Encodes every step of your data analysis as a directed acyclic graph

• Automatically detection of dependencies, which means:

  • If you change some part of your pipeline, targets will rerun (only) the parts that depend on the changed part.

Cognitive flexibility example

# _targets.R file
tar_source(files = "R") # load the R scripts in my dir
tar_option_set(packages = c("tidyverse", "marginaleffects")) # load any packages used
# define pipeline
list(
    tar_target(data_file, "data/cognitive_flexibility.csv", format = "file"),
    tar_target(data, read_csv(data_file)),
    tar_target(model, fit_model(data)),
    tar_target(marginal_effect, avg_comparisons(model, newdata = data, variables = "multilingual")),
    tar_target(model_plot, plot_predictions(model, condition = list("multilingual")))
)

tar_visnetwork()

tar_visnetwork()

Looking inside in the pipeline?

# fit_model.R
fit_model <- function(data, age_moderator = TRUE){
    if(age_moderator){
        model <- lm(cognitive_flexibility ~ multilingual * age, data = data)
    } else {
        model <- lm(cognitive_flexibility ~ multilingual, data = data)
    }
    return(model)
}

Running the pipeline

tar_make()

▶ dispatched target data_file
● completed target data_file [0.402 seconds, 17.112 kilobytes]
▶ dispatched target data
Rows: 250 Columns: 4
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (1): multilingual
dbl (3): cognitive_flexibility, age, age_c

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
● completed target data [0.114 seconds, 5.687 kilobytes]
▶ dispatched target model
● completed target model [0.003 seconds, 17.378 kilobytes]
▶ dispatched target marginal_effect
● completed target marginal_effect [0.089 seconds, 41.914 kilobytes]
▶ dispatched target model_plot
● completed target model_plot [0.036 seconds, 253.351 kilobytes]
▶ ended pipeline [0.809 seconds]

Re-inspect the pipeline

tar_visnetwork()

See what happens if we change the data

# remove anomalous observations
problematic_obs <- c(3, 15, 27)

write_csv(data[-problematic_obs,], "data/cognitive_flexibility.csv")

tar_visnetwork()

See what happens if we change the data

tar_outdated()

[1] "marginal_effect" "data"            "model_plot"      "data_file"      
[5] "model"          

See what happens if we change the model function

# R/fit_model.R, now with quadratic effect for age
fit_model <- function(data, age_moderator = TRUE){
    if(age_moderator){
        model <- lm(cognitive_flexibility ~ multilingual * poly(age, 2), data = data)
    } else {
        model <- lm(cognitive_flexibility ~ multilingual, data = data)
    }
    return(model)
}

See what happens if we change the model function

tar_visnetwork()

See what happens if we change the model function

tar_outdated()

[1] "marginal_effect" "model_plot"      "model"          

Real world example

• Big latent, multilevel timeseries model (>10,000 parameters)
• Many outputs (effect sizes, constrasts, model visualizations, data dashboard)
• Entire pipeline takes ~6 days to run

Targets at scale

Targets at scale

library(targets)
source("R/functions.R")

# Dependencies for this pipeline
tar_option_set(packages = c(
  "tidyverse",
  "rethinking",
  "rstan",
  "patchwork",
  "ggsci",
  "RColorBrewer",
  "cubelyr",
  "proxy",
  "foreach",
  "doParallel",
  "ggdist"
))

list(
  # Extract raw data
  tar_target(rdata_file,
             "data/input-to-models.RData", format = "file"),
  tar_target(footage_file,
             "data/Study1_footage_details.xlsx"),
  tar_target(infants_file,
             "data/Infant_age.xlsx"),
  tar_target(d_raw_time,
             extract_raw_data(rdata_file, footage_file, infants_file, "time")),
  
  ## Generate data for Shiny dashboard #####
  tar_target(d_dash,
             dashboard_data(d_raw_time)),

  # Data for time series model 
  tar_target(d_time,
             time_data(d_raw_time)),
  
  ## Generate data for time-averaged model comparison
  tar_target(d_avg,
             time_data_avg(d_time)),
  
  ###### Time Averaged Analyses #########################
  tar_target(model_compare,
             avg_model_comparison(d_avg)),
  
  ###### Time Series Analyses ##########################
  # Fit time series model
  tar_target(fit_time,
             fit_time_series(d_time, n_iter = 1000, n_chains = 8)),
  # Make predictions for time series
  tar_target(pred_time,
             pred_species_pid_time(fit_time, d_time)),
  
  # Make plots + export summaries of posterior predictions to .csv
  # time series plots
  tar_target(p_species_time_agent,
             plot_species_time_AOI(pred_time, "Agent")),
  tar_target(p_species_time_patient,
             plot_species_time_AOI(pred_time, "Patient")),
  tar_target(p_species_time_other,
             plot_species_time_AOI(pred_time, "Other")),
  tar_target(p_species_time_OR,
             plot_species_time_OR(pred_time)),
  tar_target(p_species_time_APO,
             plot_species_time_OR_APO(pred_time)),
  tar_target(p_species_time_AO,
             plot_species_time_OR_AO(pred_time)),
  # AOI size plots
  tar_target(p_species_AOI_size,
             plot_species_AOI_size(fit_time, d_time)),
  # AP movement diff
  tar_target(p_species_AP_md,
             plot_species_AP_md(fit_time, d_time))
  
)

Advanced use of targets

• Dynamic branching patterns, when there are too many targets to specify manually

• Great support for parallel processing using future

Conclusion / Q&A